An Inverse— Problem Approach to Designing Photonic Crystals for Cavity QED 
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Photonic band gap (PBG) materials are attractive for cavity QED experiments because they 
provide extremely small mode volumes and are monolithic, integratable structures. As such, PBG 
cavities are a promising alternative to Fabry-Perot resonators. However, the cavity requirements 
imposed by QED experiments, such as the need for high Q (low cavity damping) and small mode 
volumes, present significant design challenges for photonic band gap materials. Here, we pose the 
PBG design problem as a mathematical inversion and provide an analytical solution for a two- 
dimensional crystal. We then address a planar (2D crystal with finite thickness) structure using 
numerical techniques. 



I. INTRODUCTION 

Engineering new materials to meet specific design ob- 
jectives often begins by trial and error. New structures 
are repeatedly proposed and characterized, and the re- 
sults from each iteration are used to further refine the 
design. This process continues until an apparent opti- 
mum is achieved — that is, when the incremental modi- 
fications stop leading to improvements. 

However, from an implementation standpoint, trial 
and error is inefficient and costly, even when the process 
can be computationally simulated. Conceptually, trial 
and error provides little information about the quality of 
the optimum. This is because the design space is often 
too large to permit an exhaustive search. Therefore, it is 
common to fall back upon physical intuition (that might 
not apply to the new material) to guide the engineering 
process. Of course, this is not to say that design-by-trial 
is ineffective, only that it lacks a certain degree of rigor. 

Going beyond incremental design procedures requires 
an algorithmic, rather than intuitive, process. In many 
cases, posing the design problem as a mathematical 
inversion [|| || can provide an assessment of the resulting 
optimum. Ideally, algorithmic searches might uncover 
alternatives in the design space that physical intuition 
failed to recognize. However, such an unconstrained op- 
timum structure might prove too difficult to manufac- 
ture, in which case, the inversion optimization can be 
restricted to account for limitations in the fabrication 
capabilities. Both alternatives are beneficial. The un- 
constrained inversion provides an indication of the abso- 
lute optimal performance of the material, while the con- 
strained inversion produces the best structure that can 
actually be constructed. 

Engineering the optical properties of photonic band 
gap (PBG) structures H |, [| |6| is a process that can 
benefit from inversion techniques. Here, the objective is 



to tailor the electromagnetic modes of the crystal by ad- 
justing its spatially dependent dielectric function. Specif- 
ically, by introducing a defect into an otherwise periodic 
crystal, it is possible to produce localized electromag- 
netic fields (7). Both the spatial and temporal properties 
of these cavity modes are affected by the geometry of the 
defect. 

The ability to localize light fields has made photonic 
crystals attractive for experiments in cavity QED|| ||] 

PBG cav- 
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and the quantum information sciences) 10 
ities offer a number of advantages that make them an at- 
tractive alternative to Fabry-Perot resonators. Most no- 
tably, a small cavity mode volume is an important factor 
in achieving the strong coupling limit between trapped 
atoms and the cavity light field. Photonic crystal cav- 
ities are capable of mode volumes on the order of the 
cubic wavelength of the light 0, 12 . 
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Modern lithographic techniques should enable the in- 
tegration of PBG structures with micron-scale magnetic 
traps for neutral atoms p3| . Atom trapping experiments, 
however, pose non-trivial design challenges for photonic 
crystals. Practical concerns in cavity QED experiments, 
such as atom delivery and confinement, stron gly su ggest 
using planar photonic crystals j| ||, [3J, [3]§]l7|, []Jj, jl§] 
(two-dimensional lattices with finite thickness) rather 
than full 3D materials. However, two-dimensional PBG 
crystals only provide incomplete, quasi-3D light trap- 
ping. While well confined within the lattice plane, the 
cavity field can decay in the out-of-plane direction j|, |2(], 
^lj by coupling to radiated modes. Radiation loss should 
generally be the most significant decay mechanism in pla- 
nar photonic cavities^. Therefore, maximizing the cav- 
ity quality factor, Q = Auj/tUj, requires that this ra- 
diation loss be minimized. Simultaneously, cavity QED 
experiments require that the cavity mode function have 
high relative field strengths in vacuum regions (as op- 
posed to inside the semiconductor) that are accessible to 
trapped atoms. Otherwise, the atomic system will not 
couple strongly to the cavity field. Additionally, these 
criteria must be met without sacrificing mode volume, 
i.e., by delocalizing the defect field. 

There has recently been considerable progress to- 
ward PBG cavities that display the necessary prop- 
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erties for QED. Numerical design work performed in 
the Scherer group ||, |2| has identified planar photonic 
crystals |f8|, ^l], |23j with mode volumes on the order of 
the cubic wavelength of the light and cavity Q factors 
~10 , which is sufficient for strong coupling. However, 
with these structures, it is difficult to know if they are 
true optima. Furthermore, if they are not optimal, it is 
unclear how to further improve them, i.e., whether small 
modifications or large changes in the crystal would be 
needed. 

This paper poses the photonic crystal design process 
as an inverse problem in order to provide an algorithmic 
optimization. That is, we represent the defect dielectric 
function as resulting from the design requirements, rather 
than proposing a defect geometry and then characteriz- 
ing the crystal to observe its properties. To the best of 
our knowledge, this paper represents the first attempt 
to treat photonic band gap materials in such a man- 
ner. We demonstrate that mathematical inversion leads 
to photonic crystals structures not previously suggested 
by intuitive or trial and error design techniques, and we 
also illustrate how fabrication-imposed constraints can 
be placed on the inversion optimization. 

When solving inverse problems, there are two possible 
directions to follow. The first is to analytically treat a 
simplified model that captures the relevant properties of 
the actual problem. Analytic solutions are rarely pos- 
sible for structures of arbitrary (realistic) complexity; 
however, they provide a closed mathematical description, 
and hopefully a better understanding, of the design pro- 
cess. In Section [II we present analytic results for a pure 



(i.e. infinitely thick) two-dimensional photonic crystal. 
The second class of inversion algorithms utilize numer- 
ical methods to treat more realistic descriptions of the 
underlying physics. However, in exchange for the more 
realistic model, it can be difficult to find global extrema 
in the design space. In Section |v| , we employ numerical 
methods to treat a planar photonic crystal membrane. 



II. CAVITY DESIGN AS AN INVERSE 
PROBLEM 

The relationship between the spatially dependent di- 
electric function, ft(r)^3), and the properties of inter- 
est, such as mode volume and Q, is a composition of 
two individual components. First the electric and mag- 
netic fields are related to the reciprocal dielectric func- 
tion, rj(r) = l//t(r), through the Maxwell curl operator, 



V x 7?(r)V x Hj-(r) 



(1) 



where Hj(r) is the magnetic field for the mode, j, with 
frequency, u>j. It is often most convenient to work with 
H because the resulting Maxwell equation is Hcrmitian. 
This provides no difficulty because the electric field, 



can always be found from the magnetic field. 

The second relationship then connects the photonic 
crystal's electromagnetic modes to its characteristics. In 
some cases, the property of interest can be directly ex- 
pressed. This is true for the volume p4|, 



V, = J |^(r)| 2 dr 



(3) 



where j is the mode index and ip( r ) is defined, 

H J -(r)=H 0) ^(r) (4) 

by choosing Hp , such that tf>j is max-one normalized, 
\ip{r)\max = f@- It is also possible to directly ex- 
press the magnitude of the field at the location, r Q , of 
the trapped atom, 



h = l H j( r a = 0)| 5 



(5) 



where r a = can always be achieved by a suitable choice 
of axes. 

For other characteristics, such as the cavity Q, it may 
be too difficult, or inappropriate, to analytically express 
the property as a function of the electromagnetic modes. 
For the cavity Q, it is more convenient to work with some 
other measure, L, 



Q^L^r)] 



(6) 



that acts as a proxy for an actual calculation of the cav- 
ity Q. This measure must display the property that max- 
imizing it simultaneously maximizes the Q (this will be 
discussed in greater detail in Section III C| ) . 



A. Inversion Cost Functional 

In all these cases, the fundamental relationships that 
connect the properties of the electromagnetic modes to 
the reciprocal dielectric function, rj(r), remain implicit, 

Here, the notation, [•••], represents the fact that the 
quantities are complicated functionals of their input. 
This is because Eq. ([!]) has been buried inside them. 
Therefore, evaluating the functionals for any given crys- 
tal entails solving Maxwell's equations and then comput- 
ing the property from the resulting modes. 

Nonetheless, with these implicit functionals in hand, 
it is possible to formally state the inverse problem by 
defining a cost functional, 

J[r)(r)] = Q ro fo(r)] + z I m [v(r)] - /?vV m [r?(r)] (8) 

evaluated for the appropriate cavity mode, j = m p^| . 0i 
and (3\> are scalars that balance the relative importance of 
the various terms in the cost. Solving the inverse problem 
is accomplished by optimizing J , 



Ej(r) = iV x Hj(r)/ajje «(r) 



(2) 



maxJhfr)] 

))(r) 



(9) 
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over the possible structures (indexed by their dielectric 
function). 

Of course, Eq. (|J) is deceptively simple — all the details 
of solving the inverse problem have been relegated to the 
[■ ■ •] notation. However, there is an important advantage 
to such an abstraction. Equation (^]) provides a language 
for describing the photonic crystal structure in terms of 
the design objective. It also provides a description of the 
inversion that is independent of the particular method 
used to solve Maxwell's equations. J[r](r)] is easily gen- 
eralized to design objectives other than the cavity Q and 
mode volume, and it applies to 2D, planar, and full 3D 
materials. 

In the next two sections, Eqs. (||) and (||) are solved 
for specific examples. First, an analytic approach is used 
to treat an infinitely thick two-dimensional crystal. In 
this case, handling the [• • •] calls for the majority of the 
effort. But once this is accomplished, the optimization 
is relatively straightforward. The second example incor- 
porates numerical methods to treat a planar 2D crystal. 
Here, there is no struggle with the notation — we just 
write a computer program to compute the mode volume 
and cavity Q. However, the optimization is complicated 
by the possibility of local minima. In both cases, it is 
necessary to identify the best choice for the weighting 
parameters, Pi and /?v, in the cost functional. 



III. ANALYTICAL INVERSION 

In this section, the photonic crystal inversion is ana- 
lytically performed for a two-dimensional structure. The 
2D problem is motivated by the fact that planar and 
2D structures share many similarities. Treating the 2D 
crystal allows a detailed mathematical inversion and can 
provide insight into how to also optimize a planar crystal. 

The general inversion strategy is to solve a variational 
problem by expanding the cavity field in the bulk crys- 
tal electromagnetic modes. It is therefore necessary to 
select a bulk 2D lattice with a band gap surrounding the 
desired cavity resonance frequency (such as an hexago- 
nal array of holes with a suitable lattice constant) prior 
to the inversion. Once the electromagnetic modes of the 
bulk structure are determined, it is possible to optimize 
the cavity Q, field intensity at r a , and mode volume over 
the bulk mode expansion coefficients. This optimiza- 
tion stage does not directly involve the defect dielectric 
function — it identifies the optimal cavity field that can 
be produced using the bulk crystal modes as a basis. 
Once these optimal expansion coefficients are identified, 
the defect that produces the optimal field is extracted by 
inverting the Maxwell curl equation, Eq. (pi). 



However, it is useful to briefly review the plane wave 
expansion method in order to provide sufficient context 
for solving the inverse problem. 

The two-dimensional photonic crystal consists of a 
bulk medium with index of refraction, ni,. It is laced with 
a lattice of infinitely deep cylindrical holes with radius, 
r^, and index of refraction, n^. This lattice is repre- 
sented by the real space vectors, {R n }, which point from 
the origin to the centers of the cylinders. All of the real 
space lattice vectors lie in a plane. 

Since the inverse dielectric function, 770 (r) , for the 
bulk crystal is periodic, it is convenient to work with 
its Fourier transform, 

77o(r) = ]>> Ge 4Gr (10) 

G 

where the reciprocal lattice vectors, G, satisfy G • R n = 
2lir, I = 1,2,.... Physically, each reciprocal lattice vec- 
tor is the wave vector of a plane wave that shares the 
periodicity of the real space lattice. 

As with the dielectric function, the bulk crystal elec- 
tromagnetic modes are periodic in the lattice. In accor- 
dance with Bloch's theoremp^|, the bulk crystal electro- 
magnetic modes can also be expanded in the reciprocal 
lattice vectors, 

H„, q (r) = ^ E Vq+c^ (q+G) - r (11) 

A G 

where the mode is labelled by its wavevector, q,|36| and 
band index, n. The et\ are orthogonal polarization vec- 
tors and the ft. ra , q +G are the plane wave expansion coef- 
ficients that produce the mode. 

Calculating the bulk crystal modes is accomplished by 
solving the Maxwell equation, Eq. ([[]), using the form 
in Eq. (|TT|) . This leads to wave equations for the two 
possible polarizations, 

uj 1 

$>G-G<(q+GMq+G')V q +G< = ^rWc (12) 

G' C 

for TE modes and, 

uj 2 

E 7 ?G-G'|q + G llq + G'|v q+G , = ^f^Vq+c (13) 

G' C 

for TM modes. Since the polarizations uncouple for a 
pure two-dimensional crystals, it is possible to work with 
them independently. For the remainder of this section, 
we utilize TE modes; however, the same inversion tech- 
nique applies equally well to TM modes. 



A. Bulk Crystal Modes 



B. Defect Crystal Modes 



Methods for solving the Maxwell equations for a two- 
dimensional photonic crystal are well established M, p5|. 



The cavity mode, m, can be expanded in the bulk 
modes, H„ !q (r), using wavevectors that are confined to 
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the first Brillouin zone (see for example p6|), 

H ™«=E E ^'Hn,,(r) (14) 

of the two-dimensional lattice. This superposition is the 
reason for working with the magnetic, rather than elec- 
tric, fields. Since Eq. (Q) is Hermitian, the bulk modes 
are complete. 

The photonic crystal cavity can be described by intro- 
ducing an additional defect term, Sr](r), into the recipro- 
cal dielectric function, 

V( r ) = Vo{r) +5i](r) (15) 
whose Fourier transform is given by, 

^(r)= (16) 

A k 

However, unlike the bulk lattice, the cavity is not peri- 
odic, so the Fourier expansion must run over all wave 
vectors, k. In practice, the integral is generally approxi- 
mated by a discrete sum that is then truncated to allow 
computation. 

The coefficients, are calculated by substituting 

Eqs. (U) and @ into the Maxwell equation, Eq. (§). 
This produces the matrix eigenvalue equation, 

V D {m) a (m) = ^™ a {m) (17) 

n' ,q' 

where 

^i^n'.q' = E /l n,q+G<%q+G-q'-G'/ln',q'+G' (18) 
G,G' 

(q + G).(q / + G , )+Wq.q'^ 

Here, it can be seen that the point defect couples all of 
the bulk crystal modes. However, the Fourier coefficients, 
<5?7k, fall off quickly as the magnitude of the wavevector 
increases. 



perturbative, so there is no requirement on the relative 
magnitudes of the two functions. However, in practice, it 
can be practical to restrict the structure of the defect (for 
example, to enforce radial symmetry) in order to limit the 
number of plane waves (or equivalently, the number of 
reciprocal lattice vectors) needed for Eq. ( [j~l|) to converge. 

1. Cavity Mode Volume and Intensity 

Expressing the cavity mode volume in terms of the 
basis function coefficients is straightforward, 

Vra = EES^^'.'WI^W) ( 19 ) 
n'.q' n.q 

as is taking derivatives with respect to the coefficients, 

^ n,q 

Here, the spatial functions, il)j(r), are the quantities de- 
fined in Eq. (|) and the inner product, (Vv,q' IVVq)' de- 
notes integration over the real-space domain in which the 
mode volume is to be minimized. This domain is arbi- 
trary, provided that it is at least as large as the photonic 
crystal to be used in the cavity QED experiment. 

Similarly, the field intensity at the location of the 
trapped atom is given by, 

J ™ = E E «S<q ) H;,,q,(0)H n , q (0) (21) 
n',q' n.q 

and the necessary derivatives with respect to the expan- 
sion coefficients are also straightforward to find. When 
performing an inversion calculation, both the mode vol- 
ume and intensity functions can be further expanded in 
terms of the reciprocal lattice vector plane waves. Do- 
ing so leads to algebraic expressions in the coefficients, 
Vq+G, from Eq. (|ll|). 

2. Cavity Q Factor 



C. Cavity Mode Optimization 

To optimize the inversion cost functional, J , it is nec- 
essary to express Eq. (||) in terms of the defect crystal 
modes by utilizing the expansion, Eq. (|l4|). When per- 
forming an actual inversion calculation, this is the point 
when it would be necessary to select a bulk photonic crys- 
tal geometry, such as a hexagonal lattice. Once this has 
been done, the plane wave representations of the opti- 
mization basis functions, H„ iq (r), can be computed. 

Another important point is that the optimization is 
performed over <5?7(r), not the bulk lattice function, ?/o( r )- 
In principle, this does not restrict the optimization in 
any way. The distinction between rjo(r) and Srj(r) is not 



It is not as clear how to represent the cavity Q in 
terms of the basis functions. The two-dimensional lat- 
tice is infinitely deep and therefore does not permit any 
radiated modes. However, this does not prevent us from 
minimizing features of the two-dimensional lattice that 
would promote out-of-plane loss were the structure a pla- 
nar crystal. In other words, we wish to find the two 
dimensional structure of a planar photonic crystal that 
minimizes the out-of-plane loss by considering features 
computed for a pure 2D structure. Since the cavity Q 
cannot be directly computed, it is necessary to define 
an auxiliary measure of field decay that applies to the 
two-dimensional lattice. 

The essential requirement of the auxiliary measure is 
that optimizing it simultaneously maximizes the cavity Q 



5 



for a planar structure (where radiation loss can occur) . It 
is possible to identify such a measure by considering the 
physical nature of radiative field decay in a planar pho- 
tonic crystal. Out-of-plane loss is the result of guided 
crystal modes coupling to free-space [|| and frequency- 
wavevector pairs in free space must lie within the light 
cone. Therefore, bulk crystal modes with frequency- 
wavevector pairs that lie below the light line should not 
couple to free space because they undergo total internal 
reflection [||. 

Minimizing the contributions from bulk modes, 
H„ )q (r) that lie above the light line reduces radiative 
cavity decay. We chose to adopt the following auxiliary 
function, 



L 



/ j u n',q' u n,q 



n',q' n,q 



(22) 



where qj is the magnitude of its corresponding wavevec- 
tor, qj = \\qj\\. 



3. Analytic Optimization 

The photonic cavity design characteristics, Eqs. ( [l9f - 
|22| ), can be substituted for their respective terms in the 
cost functional, J[8rf\. Setting the derivatives of the de- 
sign properties with respect to the expansion coefficients 
equal to zero produces a linear variational problem. In 
order to insure that the mode functions remain prop- 
erly normalized, it is convenient to impose the constraint, 
1 — ^2 |o.n,q| 2 = 0, as a Lagrange multiplier. Maximiz- 
ing the resulting Lagrangian leads to a matrix eigenvalue 
problem, 



En',,' [^ L ^+PlK'^)H n ^0) 
-Pv{^n'.,q'\lpn,q)} £1^,/ = Aa„, q 



(23) 



whose eigenvectors correspond to values of the expansion 
coefficients, a^q, that satisfy Eq. (Q). The eigenvector 
corresponding to the smallest eigenvalue is the best op- 
timum. 

In practice, it is necessary to select values for the 
weighting parameters, 0j and /3y. Although appropriate 
values can be found by balancing the relative magnitudes 
of the three terms in the cost, it is better to nest Eq. 
( p4| ) within a second, numerical maximization, over the 
Pi. Although, performing the optimization calls for 
repeating the eigenvalue problem, possibly many times, 
saving the values of the variational matrix elements and 
simply rescaling them according to the particular choice 
of the (3i offers considerable computational savings. 



D. Extracting the Defect Dielectric 



With the optimal mode coefficients, known, the 

final component of the inversion is to extract the defect 



dielectric, Srj(r), from the expansion coefficients. Doing 
so involves inverting the Maxwell equations, and can be 
accomplished by substituting Eqs. ([l0|), ( pi] ) and ( |l6| ) 
into Eq. (|l|) and solving for the Sr]^. 

In simplifying the resulting expressions, it is necessary 
to make use of the orthogonality of the bulk crystal mode 
functions. It is also helpful to let the mode wavevectors, 
q, run over multiple Brillouin zones. Doing so leads to 
more manageable equations because the summations are 
no longer restricted. In the end, the proper expressions 
can be obtained by folding the equations back into the 
first Brillouin zone. The details of the derivation are 
provided in Appendix |A], and the result is a linear system 
of equations, 



E 



D 



n,q;k°'/k — «n,q 



n,q 



where the inversion matrix, T)^ m \ is given by 

^n,q;k = E E a >V,q'^rLq+G^n',q'+G-k' 
n' G,q' 

(q + G)-(q+G-k') 



(24) 



(25) 



The matrix is indexed by the bulk modes, labelled by 
(n, q), and the Fourier coefficients of the defect, k. 

An important point to make is that the cavity reso- 
nance frequency, tu m , enters into Eq. ( p4| ) as a parame- 
ter. Solving the inversion requires specifying the cavity 
frequency, which can take on any value within the bulk 
crystal band gap. It should be expected that the best in- 
plane confinement results from a cavity frequency, w m , 
deep within the band gap. However, the resulting in- 
verted defect dielectric function is different depending on 
the choice of the resonance frequency. Moreover, differ- 
ent choices of oj m might lead to defects that are easier to 
fabricate than others. This property of the inversion is 
one that likely requires further investigation. 



1. Computational Complexity 

Interpreting the computational complexity of the ma- 
trix equations in Eq. ( |25| ) is aided by considering the k 
vectors as a sum of reciprocal lattice vectors and Brillouin 
zone vectors, k = q + G. This shows that the defect is 
constructed using bulk crystal modes for all wavevectors 
that lie within the Brillouin zone. Counting vectors in 
this manner allows the dimensions of the inversion ma- 
trices to be determined. Since Eq. (|l|) is Hermitian, the 
number of bands is equal to the number of reciprocal 
lattice vectors, Nq — oo. Similarly, the number of bulk 
crystal modes needed to construct a non-periodic defect 
is N q = oo. Therefore, the dimension of the square ma- 
trix, D( m ), is a strictly countable double infinity. 

In practice, summations over reciprocal lattice vectors 
are truncated to a finite number, Nq 7 which likewise lim- 
its the number of bands for each wavevector to Nq ■ The 
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truncated bulk mode basis is no longer rigorously com- 
plete. However, in practice, the Fourier coefficients in Eq. 
( |j"o| ) fall off quickly as ||G|| increases. Therefore, truncat- 
ing the basis can be achieved without sacrificing accuracy 
provided that Fourier expansions remain sufficiently ac- 
curate. Retaining the proper dimension of D^ m \ requires 
that the number of defect Fourier wavevectors, Nk, also 
be properly chosen. The requirement that D*-™ 1 ' remain 
square requires that the number of distinct bulk modes in 
the Brillouin zone, q (recall that k = q + G) be exactly 
N G . 

The number of matrix elements in D' 1 ™* 1 scales as 
0(Nq). However, the computational complexity of the 
sums that must be evaluated when computing the ma- 
trix elements scale as 0(Nq). Therefore, the overall time 
complexity for solving the inverse problem is O(Nq), and 
can prove to be computationally demanding. In practice, 
it is beneficial to employ approximate methods for solv- 
ing the linear system of equations ]27j if the calculations 
are to be performed on a personal computer. 



E. Illustration: Symmetric Defect 

As a first demonstration, we considered an hexagonal 
2D photonic crystal with a radially symmetric defect. 
The bulk material index of refraction was chosen to be, 
n& = 3.4, and the hole radius was, r^/a = 0.3 (where 
a is the lattice constant). The rationale behind these 
parameters was based on competing factors. It has been 
shown that large lattice holes lead to out-of-plane loss 
due to scattering from the edges (|, ^2|. However, the 
frequency band gap decreases as the hole size is reduced, 
so rh must not be made too small. 30% of the lattice 
constant provides a good compromise. 



FIG. 1: The photonic crystal (A) and the cavity electric field 
(B) that result from optimizing the Q, mode volume, and peak 
cavity intensity by solving a two-dimensional inverse problem 
for a radially symmetric defect. The resulting cavity hole 
was reduced in radius to r c — 2.1/ a, and the neighboring 
holes were also reduced in size and displaced outward. The 
most significant feature of the optimized structure is that the 
index of refraction of the bulk crystal holes was increased to 
n h = 1.9. 



FIG. 2: Extracted reciprocal dielectric function, n (r) = 
l/«;(r), found by solving the inversion problem for a two- 
dimensional hexagonal lattice with a radially symmetric de- 
fect. The most distinctive feature of the mathematically op- 
timal photonic crystal is that the bulk lattice air holes have 
an index of refraction rih > 1. 

The dielectric index of refraction, = 3.4, was cho- 
sen because it is typical of the semiconductor materials 
that might be used to incorporate photonic crystals into 



atom trapping experiments (e.g., A1.3Ga.7As). The es- 
sential property of the bulk dielectric material is that it 
not absorb light around the atomic transition frequency, 
for example, 852 nm for Cs. 

For the purpose of the calculations, the photonic crys- 
tal lattice was truncated to five layers surrounding the 
center defect at r = 0. We found that this number 
of layers provided a sufficient description of the prop- 
erties of the photonic crystal without exceeding the 
computational power of a typical desktop computer. 
The coefficient optimization was performed by adjust- 
ing the parameters, /?;, in J until the best maximum 
was achieved. Parameter optimization was performed 
using a conjugate-gradient search algorithm |28| over the 
/3j. The search required solving the eigenvalue problem 
in Eq. (p3) for different scaling parameters 28 times. 



Solving the linear system of inversion equations in III D 
produced the crystal and associated cavity mode shown 
in Figure [j] Here, the 2D mode volume was ^A 2 /4. The 
location and size of the photonic crystal holes was de- 
termined from a contour plot of the reciprocal dielectric 
function, Figure ||. The contour half-way between the 
minimum and maximum values of ?y(r) was adopted in 
order to eliminate the small peak oscillations. These os- 
cillations were the result of truncating the Fourier expan- 
sions for r]o(r and Si](r). 

The photonic crystal shown in Figure |l|(A) has several 
distinctive features. First, the hole at the location of the 
cavity is reduced in radius to r c — 0.21/a. The nearest 
neighboring holes were also reduced in radius, to approx- 
imately rh = 0.26/a, as well as outwardly displaced from 
their original locations by ~ 0.15/a. Qualitative argu- 
ments for these features, which were also observed by 
Vuckovic, et al., have been proposed^. Reducing the 
size of the defect hole draws a bulk mode from the air 
band into the photonic band gap. Air band modes are 
characterized by higher intensities in regions where the 
index of refraction is lower, such as in the air hole. There- 
fore, the reduced cavity radius is most likely the result of 
the design requirement that maximizes the intensity at 
r a = 0. 

Decreasing the radii of the neighboring holes and mov- 
ing them away from the cavity reduces the intensity of 
the secondary lobes surrounding the main peak in the 
cavity mode [see Figure 0(B)]. These lobes, which would 
normally coincide with air holes in the hexagonal lattice, 
experience an atypically high index of refraction. The 
lobe intensities are suppressed because the defect mode, 
which was pulled from the air band, is low-index seeking. 
As a result, the cavity mode displays better localization. 
Therefore, the displaced neighboring holes likely result 
from the mode volume minimization design requirement. 

By far, the most significant feature of the optimized 
photonic crystal in Figure [l] is that the index of refrac- 
tion of the bulk holes is increased to ~ 1.9. The opti- 
mal structure is a bulk photonic crystal whose holes are 
made of a material other than air. Understanding this 
result requires analyzing how the cavity mode is con- 
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FIG. 3: Dispersion relationship and bulk mode contributions 
to the optimized cavity mode for a radially symmetric defect. 
The upper shaded region denotes the free space light line. The 
band structure (solid lines) of the optimized photonic crystal 
is shown and the points indicate bulk modes that contribute 
more than 1% to the cavity mode. The dotted lines represent 
the band structure of a photonic crystal with air holes. 



structed from the bulk mode basis functions. The solid 
lines in Figure [| show the dispersion relationship of the 
optimized crystal (with non-air holes). For reference, the 
dispersion relationship for a corresponding crystal with 
air holes {rih = 1.0) is denoted by the broken lines. The 
points (circles) represent bulk modes, identified by their 
band index, wavevector and frequency, that contribute 

more than 1%, |al"q| 2 > 10~ 2 , to the optimized cavity 
mode. 

As can be seen in Figure the crystal with rih = 1.9 
contains more bands that lie below the light line. Con- 
sequently, more modes, particularly with larger wavevec- 
tors and frequencies, contribute to the cavity mode with- 
out sacrificing Q (since these modes now lie below the 
light line). Most likely, the ability to include contribu- 
tions from more bulk crystal modes allows an increase in 
the cavity Q while simultaneously decreasing the mode 
volume. 

The rationale here is that leaky modes must be ex- 
cluded by the Q maximization. However, constructing a 
cavity mode that is well localized, i.e., not periodic, re- 
quires a large number of Fourier components, including 
high frequency modes. Drawing bands out of the light 
cone provides access to more basis functions and this al- 
lows a localized cavity field to be constructed without 
resorting to leaky modes. Still, the optimization did in- 
corporate several bulk modes from above the light line, 
particularly from the J — T line. 

We do not attempt to comment on how photonic crys- 
tals with non-air holes might be fabricated. However, 
it is important to consider the major effects of adjusting 
the hole dielectric. Increasing the hole index of refraction 
decreases the band gap. However, because the index con- 
trast between the bulk material and the holes is smaller, 
it should be possible to increase the hole radius without 
suffering as much vertical scattering from the edges. Of 
course, increasing the hole size pushes bands back into 
the light cone and a balance must be found. 



F. Illustration: Asymmetric Defect 

As a second demonstration, we relaxed the radial sym- 
metry requirement, and considered an arbitrary defect 
in a two-dimensional hexagonal latti ce. T he same pho- 



tonic crystal parameters from Section III E were adopted: 
a bulk index of refraction, nj — 3.4, and hole radius, 
fh/o- = 0.3. Again, the mode expansion was constructed 
using five photonic crystal layers surrounding the defect 



center to provide a sufficient bulk mode basis expansion 
without incurring excessive computational expense. 



FIG. 4: The photonic crystal (A) and the cavity mode (B) 
that result from optimizing the Q, mode volume, and cavity 
intensity by solving a two-dimensional inverse problem (radial 
symmetry not imposed). The cavity hole was reduced in ra- 
dius to r c = 2.18/a, and the neighboring holes were elongated 
in the vertical direction. The most significant feature of the 
optimized structure is that the index of refraction of the bulk 
crystal holes was increased to n/, = 1.75. 



FIG. 5: Dispersion relationship and bulk mode contributions 
to the optimized cavity mode for a radially symmetric defect. 
The upper shaded region denotes the free space light line. The 
band structure (solid lines) of the optimized photonic crystal 
is shown and the points indicate bulk modes that contribute 
more than 1% to the cavity mode. The broken lines represent 
the band structure of a photonic crystal with air holes. 

Solving the linear system of inversion equations 
was performed by nesting the optimization within a 
conjugate-gradient search for the best Pi weighting pa- 
rameters. The photonic crystal was again constructed 
by taking the contour level half way between the mini- 
mum and maximum of rj(r) to eliminate the numerical 
effects of truncating the Fourier expansions. 

The resulting photonic crystal and cavity mode is de- 
picted in Figure Here, the volume of the 2D mode 
was reduced to ~"/10. The center photonic crystal hole 
was reduced in radius to r c = 0.22/a and remained circu- 
lar. However, the nearest neighbor holes were deformed 
mainly in the y direction. The minor axes, along x, of 
the four holes above and below the cavity were reduced 
to 77, ~ 0.26/a. Their corresponding major axes were 
simultaneously increased to 0.39/a, resulting in the el- 
liptical holes seen in Figure ||(A). These four neighbor- 
ing holes were also radially displaced from their normal 
hexagonal lattice sites by 0.05/a, smaller than what was 
observed for the symmetric defect. 

There was also a deformation of holes lying along the 
i-axis. These were elongated in the y direction; how- 
ever, the degree of eccentricity was not constant. The 
minor axes of the two horizontally neighboring holes de- 
creased to 0.26/a, while their major axes increased to 
0.31/a. The remaining horizonal holes were also verti- 
cally stretched, with major axes given by .36/a, .31/a 
and .30/a for photonic crystal layers 3, 4 and 5, respec- 
tively. Slight deformations along the x direction were 
observed in the y-axis holes. However, their deforma- 
tion was small, with major axes of 0.31/a along the x 
direction, and minor axes of 0.29/a. No consistent defor- 
mation of hole size was observed for any other direction. 

The qualitative arguments suggested by Vuckovic, et 
ai.Sjj also apply to this structure. Replacing air hole 
sites with higher index material suppresses the amplitude 
of air-band originated modes. But more importantly, the 
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FIG. 6: Schematic of the planar photonic crystal optimized 
using numerical inversion methods, a is the hexagonal lattice 
constant, d is the thickness of the slab and ru is the radius of 
the bulk crystal holes. The design parameters were obtained 
by computing the crystal electromagnetic fields produced by 
plane wave illumination. The cavity Q was obtained by com- 
puting the reflection coefficient, R(u>) as a function of incident 
wave frequency. 



partial elongation of holes along the x axis bears signif- 
icant resemblance to the "fractional edge delocations" 
suggested by the Scherer group |8|, ^jj. 

As with the radially symmetric defect, the index of re- 
fraction of the bulk holes increased, but only to ~ 1.8. 
A similar argument based on pulling higher bands out of 
the light cone again applies. However, it is not as easy to 
provide an argument for the elongated holes surrounding 
the defect and along the x and y axes. Nonetheless, it can 
be seen from Figure |J(B) that better mode localization 
was achieved by allowing an asymmetric cavity defect. 

Figure [| shows the band diagram for the optimized 
photonic crystal. The dispersion relationship for bulk 
holes with = 1.8 is depicted by the solid lines, and 
the air hole bands are given by the broken lines. The cir- 
cles represent bulk modes that contribute more than 1%, 
o-n^q > 10~ 2 , to the optimized cavity mode. The same 
argument used to explain the symmetric defect results 
can be made here. The crystal with — 1.8 contains 
more bands below the light line. 

A distinctive feature of the optimal asymmetric mode 
is that it contains fewer contributions from leaky modes 
(inside the light cone) than for the symmetric case. This 
can be seen by comparing Figures || and [|. Exactly how 
the larger Q was achieved for the asymmetric cavity is 
not clear. However, it is not surprising that constraining 
the optimization (for example by imposing a symmetry 
restriction) reduces the ability of the inversion to satisfy 
the design requirements. 

In these two-dimensional examples, solving the design 
problem by inversion resulted in structures that differ 
from the previous, trial and error designs that have been 
suggested. These new photonic crystal structures, how- 
ever, might prove difficult to fabricate. As such, they 
provide an indication of the performance of an ideal, per- 
haps experimentally impractical, cavity. 



IV. NUMERICAL INVERSION RESULTS 

As a final example of inverting a photonic crystal de- 
fect, we considered a hexagonal structure with finite 
depth to enable a direct calculation of the cavity Q. 
Additionally, we constrained the index of refraction of 
the crystal holes to nj, = 1. In this demonstration, the 
aim was to identify an optimal planar photonic structure 
without requiring a more complicated fabrication. 

For the planar structure, it was necessary to aban- 



don an analytical solution and perform the inversion op- 
timization numerically. In some respects, the numeri- 
cal design problem was much simpler than the analytic 
case — it was possible to maximize the inversion cost 
functional, J ', without the need for a specific mathe- 
matical analysis. Instead, the spatial dependence of the 
dielectric function was parameterized. Then the inverse 
problem was solved by optimizing t7[?7(r)] over the pa- 
rameter space using a genetic algorithm. 

However, numerical methods introduce several new 
challenges. For example, since the inversion cost func- 
tional involved a large number of parameters, the result- 
ing multi-dimensional optimization was complicated by 
the existence of local minima. Additionally, numerically 
integrating the Maxwell equations is computationally ex- 
pensive. For instance, a single mode volume and cavity 
Q calculation can require several minutes of computer 
time. 



A. Planar Photonic Crystal 

We considered an eight layer, two-dimensional hexag- 
onal lattice slab, as depicted in F igure |[ Based on the 
previous arguments (c./., Section HIE), the hole radius 
was chosen to be r/j/a = 0.3 with a bulk material index 
of refraction, = 3.4. The slab thickness was taken to 
be, d/a — 3/4 in order to prevent multimode behavior [Q. 
Parameterizing the dielectric function was accomplished 
by treating the positions and major and minor axes of 
the lattice holes as optimization variables. Specifically, 
the cavity hole, the holes surrounding the defect, and 
the i-axis holes were considered. In all, there were 52 
optimization variables, 4 for 13 different lattice sites. 

The Maxwell equations were solved by employing the 
transfer matrix method (TMM) of Pendry and Bell[p9[. 
This technique provided a means for computing the sta- 
tionary states of the electromagnetic field for a crystal 
illuminated on one edge by a plane wave (refer to Fig- 
ure H). The spatial dependence of the electromagnetic 
fields, as well as the crystal transmission and reflection 
coefficients were calculated using a first order finite dif- 
ference solution to Maxwell's equations. The integration 
mesh used to compute the wavefields extended five lay- 
ers above and below the surface of the slab and absorbing 
boundary conditions were imposed at the top and bottom 
of the integration cell. 

The mode volume was computed by integrating the 
square magnitude of the electric field over the interior of 
the photonic crystal. The volume integration was per- 
formed using a three dimensional Simpson's rule quadra- 
ture, and the result was scaled by the field maximum, ac- 



FIG. 7: Optimization progress as a function of genetic algo- 
rithm generation number for the numerical inversion involving 
a planar photonic crystal. 
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cording to Eq. (||). The cavity Q was directly computed 
by scanning the reflection coefficient over the frequency 
band gap. The full width at half maximum (FWHM) of 
the reflection lineshape was then used to determine the 
cavity Q. 

The inversion cost functional was optimized by em- 
ploying a genetic algorithm p0[ (GA) to maximize 
J[q(Y)]. A GA was chosen because, although not terri- 
bly efficient, genetic algorithms provide good exploration 
to exploitation in multi-dimensional searches. As such, 
they are generally successful at avoiding local minima on 
a complex optimization surface. However, good explo- 
ration does not come without cost. 27 GA generations 
(iterations of the algorithm) were required to optimize 
J , as can be seen Figure 0. With a population size of 
10, a steady-state propagation routine, a mutation rate 
of 15% and a crossover rate of 85%, the full optimization 
required 24.5 hours of CPU time on an Intel PHI 1.2 GHz 
desktop machine. 

The resulting photonic crystal and the optimized cav- 
ity mode are shown in Figure |^. As can be seen, a sim- 
ilar structure to the two-dimensional optimization was 
obtained. The cavity defect hole was reduced in radius 
to 0.21 1 a and the holes surrounding the defect were also 
decreased in size. The holes along the i-axis were also 
elongated in the z direction, which again demonstrated 
remarkable similarity to the fractional edge delocations 
suggested by the Scherer group (§1, p2| . 

The mode volume of the field in Figure ^ was found to 
be V m ~ A 3 /3 and the cavity Q was ~1.1 x 10 5 . These 
quantities surpass those of the best PBG cavity designs 
to date. 



V. CONCLUSION 

We demonstrated that inversion methods provide a 
powerful technique for designing photonic crystals for 
cavity QED experiments. Both an analytical solution 
for a two-dimensional crystal and numerical results for a 
planar (finite thickness) crystal were presented. In both 
cases, the design objectives, namely high cavity Q and 
small mode volume, were achieved. 

For two-dimensional crystals, it was shown that contri- 
butions to the cavity mode from leaky bulk modes could 
be minimized without delocalizing the field. Here, in- 
version produced design alternatives not previously sug- 
gested by trial and error, or intuitive design. The optimal 
cavity defect for both the radially symmetric and asym- 
metric defect optimizations called for a lattice of holes 
with an index of refraction greater than one. An ex- 
planation for this inversion result, i.e., bulk crystal holes 



composed of a material other than air, was suggested. In- 
creasing the index of refraction of the lattice holes pulls 
more of the band diagram out of the light cone. Conse- 
quently, less of the higher frequency bulk modes are leaky. 
It was argued that a better localized cavity field can be 
constructed by incorporating these higher frequency ba- 
sis functions. 

The inverse problem design approach was also applied 
to a planar crystal in order to treat a more realistic struc- 
ture. However, this required that numerical methods for 
integrating the Maxwell equations as well as optimizing 
the inversion cost function be adopted. It was demon- 
strated that a crystal with low mode volume and high 
Q could be achieved. In order to minimize the effects of 
local minima in the inversion optimization, a genetic al- 
gorithm was adopted, despite its computational expense. 

For the planar photonic crystal inversion, the dielec- 
tric constant of the bulk lattice holes was constrained 
to tij, = 1. This provided an example of restricting the 
inversion optimization because of fabrication considera- 
tions. Although it may be possible to manufacture the 
inversion designs from Section III, i.e., bulk lattice holes 



with n/i > 1, it is also important to identify the best 
possible structure which can fabricated using currently 
available techniques. 

An aspect of the design process that was not consid- 
ered in this paper was robustness. Ideally, the photonic 
crystal would be insensitive to slight variations in its 
structure, as might result from small fabrication errors, 
temperature fluctuations, and etc. Including robust- 
ness into design optimization problems would add ad- 
ditional constraints into the optimization process, how- 
ever, algorithmic approaches to finding robust optima are 
known |pl|, |32| . In addition to investigating robustness, 
we plan to explore the use of optimization methods to de- 
sign PBG structures with convenient geometric features, 
such as enlarged holes for atom trapping. 

It was shown that the mathematical analysis tools uti- 
lized here simultaneously optimize multiple, complimen- 
tary design requirements. Achieving similar outcomes via 
trial and error methods is generally a formidable task. 
Therefore, an algorithmic approach to design problems 
with multiple objectives, such as with photonic band gap 
materials, will likely only be possible via an algorithmic 
approach. 
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APPENDIX A: 2D INVERSION EQUATIONS 

In this appendix we describe in greater detail how to 
extract the defect dielectric function from the optimized 



expansion coefficients (refer to Section HID). The gen- 



eral procedure involves substituting the normal mode ex- 
pansion, Eq. (fill), and the Fourier expansions, Eqs. jl0| ) 
and ([l6|) , into the Maxwell equation, Eq. ([!]) . The result- 
ing inversion equation, Eq. (|24|), results from employing 
the orthogonality of the bulk modes. 

In doing the algebra, it is more convenient to work a 
bulk mode expansion that includes wave vectors from all 
of k-space, i.e., by incorporating multiple Brillouin zones, 



Hn 



(Al) 



+V x (fy(r)V x H m (r) = ^H m (r) 
left multiplying by, H„». q «, and integrating leads to, 

* E «S%^^ H n"^(r)H n)k (r)dr = 

n, k 

iEE E a S$mcK>>M>>+G'' h n,*+G (A3) 

k n,kG,G" 

(k" + G") • (k" + G" - k') 

J^ e i(k+G+k'-k"-G")-r dr 

where, again, the integrations are over the N unit cells. 
Applying the orthogonality of the bulk modes leads to a 
set of equations that still run over all wavevectors, 



where N is the number of Brillouin zones (reciprocal lat- 
tice vectors). The orthogonality relationships are now 
give by, 

f H*, ik ,(r)H„ )k (r)dr = 5 n , n ^^ !k+G (A2) 

where the integration is over the N associated Wigner- 
Seitz cells. 

Substituting both the bulk crystal and defect dielectric 
Fourier expansions into the Maxwell curl equation, 

V x 7j (r)V x H m (r) 



lEEE 4 m k ) ( %^,k"+G" Vc"-k< (A4) 

n,k k' G" 

( k " + g") . ( k " + g" - 10 = ^«K k » "-" k ;""" 



Then, the summations can be folded back into the first 
Brillouin zone by using the identity, 

E«iS = ( A5 ) 

n,k n,q 

and the indices can be renamed to produce Eq. (E4h . 
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